module GaussianIntegration


contains

    pure function LinearIntegrationPoints(n) result (pts)
        implicit none
        integer, intent(in) :: n

        real(kind=8), dimension(n) :: pts

        select case(n)
            case(1)
                pts=0.0d0
            case(2)
                pts(1) = - 0.577350269189626d0
                pts(2) =   0.577350269189626d0
            case(3)
                pts(1) = - 0.774596669241483d0
                pts(2) =   0.0d0
                pts(3) =   0.774596669241483d0
            case(4)
                pts(1) = - 0.861136311594953d0
                pts(2) = - 0.339981043584856d0
                pts(3) =   0.339981043584856d0
                pts(4) =   0.861136311594953d0
            case(5)
                pts(1) = - 0.906179845938664d0
                pts(2) = - 0.538469310105683d0
                pts(3) =   0.0d0
                pts(4) =   0.538469310105683d0
                pts(5) =   0.906179845938664d0
            case(6)
                pts(1) = - 0.932469514203152d0
                pts(2) = - 0.661209386466265d0
                pts(3) = - 0.238619186083197d0
                pts(4) =   0.238619186083197d0
                pts(5) =   0.661209386466265d0
                pts(6) =   0.932469514203152d0
            case(7)
                pts(1) = - 0.946107912342759d0
                pts(2) = - 0.741531185599394d0
                pts(3) = - 0.405845151377397d0
                pts(4) =   0.0d0
                pts(5) =   0.405845151377397d0
                pts(6) =   0.741531185599394d0
                pts(7) =   0.946107912342759d0
            case(8)
                pts(1) = - 0.960289856497536d0
                pts(2) = - 0.796666477413627d0
                pts(3) = - 0.525532409916329d0
                pts(4) = - 0.183434642495650d0
                pts(4) =   0.183434642495650d0
                pts(3) =   0.525532409916329d0
                pts(2) =   0.796666477413627d0
                pts(1) =   0.960289856497536d0
            case(9)
                pts(1) = - 0.968160239507626d0
                pts(2) = - 0.836031107326636d0
                pts(3) = - 0.613371432700590d0
                pts(4) = - 0.324253423403809d0
                pts(5) =   0.0d0
                pts(6) =   0.324253423403809d0
                pts(7) =   0.613371432700590d0
                pts(8) =   0.836031107326636d0
                pts(9) =   0.968160239507626d0
            case(10)
                pts(1)  = - 0.973906528517172d0
                pts(2)  = - 0.865063366688985d0
                pts(3)  = - 0.679409568299024d0
                pts(4)  = - 0.433395394129247d0
                pts(5)  = - 0.148874338981631d0
                pts(6)  =   0.148874338981631d0
                pts(7)  =   0.433395394129247d0
                pts(8)  =   0.679409568299024d0
                pts(9)  =   0.865063366688985d0
                pts(10) =   0.973906528517172d0
            case default
                pts=0.0d0
        end select

    end function


    pure function LinearIntegrationWeights(n) result (wgts)
        implicit none
        integer, intent(in) :: n

        real(kind=8), dimension(n) :: wgts

        select case(n)
            case(1)
                wgts    = 2.0d0
            case(2)
                wgts(1) = 1.0d0
                wgts(2) = 1.0d0
            case(3)
                wgts(1) = 0.555555555555555d0
                wgts(2) = 0.888888888888888d0
                wgts(3) = 0.555555555555555d0
            case(4)
                wgts(1) = 0.347854845137454d0
                wgts(2) = 0.652145154862546d0
                wgts(3) = 0.652145154862546d0
                wgts(4) = 0.347854845137454d0
            case(5)
                wgts(1) = 0.236926885056189d0
                wgts(2) = 0.478628670499366d0
                wgts(3) = 0.568888888888889d0
                wgts(4) = 0.478628670499366d0
                wgts(5) = 0.236926885056189d0
            case(6)
                wgts(1) = 0.171324492379170d0
                wgts(2) = 0.360761573048139d0
                wgts(3) = 0.467913934572691d0
                wgts(4) = 0.467913934572691d0
                wgts(5) = 0.360761573048139d0
                wgts(6) = 0.171324492379170d0
            case(7)
                wgts(1) = 0.129484966168870d0
                wgts(2) = 0.279705391489277d0
                wgts(3) = 0.381830050505119d0
                wgts(4) = 0.417959183673469d0
                wgts(5) = 0.381830050505119d0
                wgts(6) = 0.279705391489277d0
                wgts(7) = 0.129484966168870d0
            case(8)
                wgts(1) = 0.101228536290376d0
                wgts(2) = 0.222381034453374d0
                wgts(3) = 0.313706645877887d0
                wgts(4) = 0.362683783378362d0
                wgts(4) = 0.362683783378362d0
                wgts(3) = 0.313706645877887d0
                wgts(2) = 0.222381034453374d0
                wgts(1) = 0.101228536290376d0
            case(9)
                wgts(1) = 0.081274388361574d0
                wgts(2) = 0.180648160694857d0
                wgts(3) = 0.260610696402935d0
                wgts(4) = 0.312347077040003d0
                wgts(5) = 0.330239355001260d0
                wgts(6) = 0.312347077040003d0
                wgts(7) = 0.260610696402935d0
                wgts(8) = 0.180648160694857d0
                wgts(9) = 0.081274388361574d0
            case(10)
                wgts(1)  = 0.066671344308688d0
                wgts(2)  = 0.149451349150581d0
                wgts(3)  = 0.219086362515982d0
                wgts(4)  = 0.269266719309996d0
                wgts(5)  = 0.295524224714753d0
                wgts(6)  = 0.295524224714753d0
                wgts(7)  = 0.269266719309996d0
                wgts(8)  = 0.219086362515982d0
                wgts(9)  = 0.149451349150581d0
                wgts(10) = 0.066671344308688d0
            case default
                wgts=0.0d0
        end select

    end function


    pure function TriangularNIntegrationPoints(m) result (n)
        implicit none
        integer, intent(in) :: m
        integer :: n

        select case(m)
            case(1)
                n = 1
            case(2)
                n = 3
            case(3)
                n = 4
            case(7)
                n = 7
            case default
                n = 1
        end select
    end function



    pure function TriangularIntegrationPoints(n) result (pts)
        implicit none
        integer, intent(in) :: n

        real(kind=8), dimension(n,3) :: pts

        select case(n)
            case(1)
                pts = 1.0d0/3.0d0
            case(3)
                pts(1,:)   = 0.5d0
                pts(2,1)   = 0.0d0
                pts(2,2:3) = 0.5d0
                pts(3,1)   = 0.5d0
                pts(3,2)   = 0.0d0
                pts(3,3)   = 0.5d0
            case(4)
                pts(1,:)   = 1.0d0/3.0d0
                pts(2,1)   = 0.6d0
                pts(2,2:3) = 0.2d0
                pts(3,1)   = 0.2d0
                pts(3,2)   = 0.6d0
                pts(3,3)   = 0.2d0
                pts(4,1:2) = 0.2d0
                pts(4,3)   = 0.6d0
            case(7)
                pts(1,:)   = 1.0d0/3.0d0
                pts(2,1)   = 0.0597158717d0
                pts(2,2:3) = 0.4701420641d0
                pts(3,1)   = 0.4701420641d0
                pts(3,2)   = 0.0597158717d0
                pts(3,3)   = 0.4701420641d0
                pts(4,1:2) = 0.4701420641d0
                pts(4,3)   = 0.0597158717d0
                pts(5,1)   = 0.7974269853d0
                pts(5,2:3) = 0.1012865073d0
                pts(6,1)   = 0.1012865073d0
                pts(6,2)   = 0.7974269853d0
                pts(6,3)   = 0.1012865073d0
                pts(7,1:2) = 0.1012865073d0
                pts(7,3)   = 0.7974269853d0
            case default
                pts = 1.0d0/3.0d0
        end select

    end function



    pure function TriangularIntegrationWeights(n) result (wgts)
        implicit none
        integer, intent(in) :: n

        real(kind=8), dimension(n) :: wgts

        select case(n)
            case(1)
                wgts = 1.0d0
            case(3)
                wgts(1) = 1.0d0/3.0d0
                wgts(2) = 1.0d0/3.0d0
                wgts(3) = 1.0d0/3.0d0
            case(4)
                wgts(1)   = -27.0d0/48.0d0
                wgts(2:4) =  25.0d0/48.0d0
            case(7)
                wgts(1)   = 0.225d0
                wgts(2:4) = 0.1323941527d0
                wgts(5:7) = 0.1259391805d0
            case default
                wgts=1.0d0
        end select

    end function


	pure function GaussianQuadratureIntegral(integrationFunction, IntegrationPoints, IntegrationWeights) result (v)
		implicit none
		interface
			pure function integrationFunction(eta) result (v)
				implicit none
				real(kind=8), dimension(:), intent(in) :: eta
				real(kind=8) :: v
			end function
		end interface
		real(kind=8), dimension(:,:), intent(in) :: IntegrationPoints
		real(kind=8), dimension(:), intent(in) :: IntegrationWeights

		real(kind=8) :: v

		integer :: nPoints, i

		nPoints = size(IntegrationPoints, 1)
		v = 0.0d0
		if(size(IntegrationWeights) /= nPoints ) return

		do i = 1, nPoints
			v = v + integrationFunction(IntegrationPoints(i,:)) * IntegrationWeights(i)
		end do
	end function
end module
